#######################################
#Quasi-Poisson Models with Simulations
#######################################

rm(list = ls())
setwd() #Set Working Directory

# Packages
pacman::p_load(rtweet,
               ggplot2,
               purrrlyr,
               dplyr,
               readxl,
               writexl,
               broom,
               tidyverse,
               dotwhisker,
               lubridate,
               stargazer,
               ggpubr, 
               jtools,
               rpart,
               caret,
               reshape2,
               ipred,
               RColorBrewer,
               rpart.plot,
               MASS,
               randomForest,
               dplyr,
               car,
               vtable,
               effects,
               interactions,
               doParallel,
               foreach
               
)

set.seed(0611)
#Import data
df <- readRDS("df.RDS")

#QP Models
##########
##########

# Creating wartime variable 
df %>%
  mutate(War = created_at > "2022-02-24") -> df

# coercing for formatting 
df$War <- ifelse(df$War == "TRUE", 1, 0)
df$user_id <- as.factor(df$user_id)   
df$conspir.po <- as.factor(df$conspir.po)
df$conspir.collusion <- as.factor(df$conspir.collusion)
df$conspir.rus <- as.factor(df$conspir.rus)
df$conspir.tusk <- as.factor(df$conspir.tusk)
df$month <- as.factor(df$month)

#########
#########
# Models
#########
#########

##########
#Likes
##########
m.conspir.po.fav <- glm(favorite_count ~ conspir.po*War*officials + log_follow + 
                          log_friends + verified, 
                        family = "quasipoisson",
                        data = df)
#same but for stargazer
m1 <- glm(favorite_count ~ conspir.po*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)

m.conspir.tusk.fav <- glm(favorite_count ~ conspir.tusk*War*officials + log_follow + 
                            log_friends + verified + month, 
                          family = "quasipoisson",
                          data = df)
#same but stargazer 
m2 <- glm(favorite_count ~ conspir.tusk*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)

m.conspir.rus.fav <- glm(favorite_count ~ conspir.rus*War*officials + log_follow + 
                           log_friends + verified + month, 
                         family = "quasipoisson",
                         data = df)
#same but stargazer
m3 <- glm(favorite_count ~ conspir.rus*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)

m.conspir.collusion.fav <- glm(favorite_count ~ conspir.collusion*War*officials + log_follow + 
                                 log_friends + verified + month, 
                               family = "quasipoisson",
                               data = df)
#same but for stargazer
m4 <- glm(favorite_count ~ conspir.collusion*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)

stargazer(m1, m2, m4,
          #m.conspir.po.rt, m.conspir.tusk.rt, m.conspir.collusion.rt, 
          type = "latex",
          header = FALSE,
          align = TRUE,
          label = "tab:triple.inter.qp",
          title = "Tweet Characteristics and Frequency of Likes- Quasi-Poisson Models",
          dep.var.labels = c("Number of Likes"
                             # "Number of Retweets"
          ),
          model.names = FALSE,
          omit = c("month"), #so it's prettier
          covariate.labels = c("PO Consp.",
                               "Tusk Consp.",
                               "Russia and PO Consp.",
                               "War",
                               "PiS Officials",
                               "Log(N Followers)",
                               "Log(N Friends)",
                               "Verified",
                               "PO Consp.*War",
                               "PO Consp.*PiS Official",
                               "Tusk Consp.*War",
                               "Tusk Consp.*PiS Official",
                               "Collusion Consp.*War",
                               "Collusion Consp.*PiS Official",
                               "War*PiS Official",
                               "PO Consp.*War*Official",
                               "Tusk Consp.*War*Official",
                               "Collusion Consp.*War*Official",
                               "Constant"),
          notes.append = FALSE,
          notes = c("Standard Errors in parenthesis.", 
                    "*p $<$ 0.1; **p $<$ .05, ***p $<$ .01 (two-tailed tests)."))
##########
#Retweets
##########

m.conspir.po.rt <- glm(retweet_count ~ conspir.po*War*officials + log_follow + 
                         log_friends + verified + month, 
                       family = "quasipoisson",
                       data = df)
m1 <- glm(retweet_count ~ conspir.po*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)
m.conspir.tusk.rt <- glm(retweet_count ~ conspir.tusk*War*officials + log_follow + 
                           log_friends + verified + month, 
                         family = "quasipoisson",
                         data = df)
m2 <- glm(retweet_count ~ conspir.tusk*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)
m.conspir.rus.rt <- glm(retweet_count ~ conspir.rus*War*officials + log_follow + 
                          log_friends + verified + month, 
                        family = "quasipoisson",
                        data = df)
m3 <- glm(retweet_count ~ conspir.rus*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)
m.conspir.collusion.rt<- glm(retweet_count ~ conspir.collusion*War*officials + log_follow + 
                               log_friends + verified + month, 
                             family = "quasipoisson",
                             data = df)
m4 <- glm(retweet_count ~ conspir.collusion*War*officials + log_follow + 
            log_friends + verified + month, 
          family = "quasipoisson",
          data = df)

stargazer(m1, m2, m4,
  type = "latex",
  header = FALSE,
  align = TRUE,
  label = "tab:triple.inter.qp",
  title = "Tweet Characteristics and Frequency of Retweets- Quasi-Poisson Models",
  dep.var.labels = c(#"Number of Likes",
    "Number of Retweets"
  ),
  model.names = FALSE,
  omit = c("month"), #so it's prettier
  covariate.labels = c("PO Consp.",
                       "Tusk Consp.",
                       "Russia and PO Consp.",
                       "War",
                       "PiS Officials",
                       "Log(N Followers)",
                       "Log(N Friends)",
                       "Verified",
                       "PO Consp.*War",
                       "PO Consp.*PiS Official",
                       "Tusk Consp.*War",
                       "Tusk Consp.*PiS Official",
                       "Collusion Consp.*War",
                       "Collusion Consp.*PiS Official",
                       "War*PiS Official",
                       "PO Consp.*War*Official",
                       "Tusk Consp.*War*Official",
                       "Collusion Consp.*War*Official",
                       "Constant"),
  notes.append = FALSE,
  notes = c("Standard Errors in parenthesis.", 
            "*p $<$ 0.1; **p $<$ .05, ***p $<$ .01 (two-tailed tests)."))

###############
#Bootstrapping 
##############
#Create function for bootstrapping
boot_pi <- function(model, pdata, n, p) {
  odata <- model$data
  lp <- (1 - p) / 2
  up <- 1 - lp
  set.seed(2016)
  seeds <- round(runif(n, 1, 1000), 0)
  boot_y <- foreach(i = 1:n, .combine = rbind) %dopar% {
    set.seed(seeds[i])
    bdata <- odata[sample(seq(nrow(odata)), size = nrow(odata), replace = TRUE), ]
    bpred <- predict(update(model, data = bdata), type = "response", newdata = pdata)
    rpois(length(bpred), lambda = bpred)
  }
  boot_ci <- t(apply(boot_y, 2, quantile, c(lp, up)))
  return(data.frame(pred = predict(model, newdata = pdata, type = "response"), lower = boot_ci[, 1], upper = boot_ci[, 2]))
}
###
#PO
###
sim.df <- data.frame(conspir.po = sample(rep(c(0,1,0,1), length.out = 250)),
                     War = sample(rep(c(0,1,0,1), length.out = 250)),
                     officials = sample(rep(c(0,1,0,1), length.out = 250)),
                     log_follow = sample(seq(min(df$log_follow, na.rm = T), 
                                             max(df$log_follow, na.rm = T),
                                             length.out = 500)),
                     log_friends = sample(seq(min(df$log_friends, na.rm = T), 
                                              max(df$log_friends, na.rm = T),
                                              length.out = 500)),
                     verified = sample(c(rep(TRUE, 33), rep(FALSE, 467))),
                     month = c(rep(1:3,length.out = 42), 
                               rep(5:12, length.out = 112),
                               rep(4, length.out = 346)))

sim.df$conspir.po <- as.factor(sim.df$conspir.po)
sim.df$month <- as.factor(sim.df$month) 

preds <-  predict.glm(m.conspir.po.fav, sim.df, type = "response", se.fit = F)
preds2 <- cbind(sim.df, boot_pi(m.conspir.po.fav, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.po, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))


#Calculate differences

po.fav.diff <- as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                            preds3$upper2[5] - preds3$upper2[1],
                                            preds3$lower2[5] - preds3$lower2[1])), 
                                     0, #War =0
                                     0), #officials = 0
                                   c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                            preds3$upper2[6] - preds3$upper2[2],
                                            preds3$lower2[6] - preds3$lower2[2])),
                                     0, #War =0
                                     1),  #officials = 1
                                   c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                            preds3$upper2[7] - preds3$upper2[3],
                                            preds3$lower2[7] - preds3$lower2[3])),
                                     1, #War =1
                                     0), #officials = 0
                                   c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                            preds3$lower2[8] - preds3$lower2[4],
                                            preds3$upper2[8] - preds3$upper2[4])),
                                     1, #War =1
                                     1) #officials = 1
                                   
))

colnames(po.fav.diff) <- c("LCI", "estimate", "UCI", "War", "PiS")


#PO & Retweets

preds2 <- cbind(sim.df, boot_pi(m.conspir.po.rt, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.po, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))

#Calculate differences
po.rt.diff <- as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                           preds3$upper2[5] - preds3$upper2[1],
                                           preds3$lower2[5] - preds3$lower2[1])), 
                                    0, #War =0
                                    0), #officials = 0
                                  c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                           preds3$upper2[6] - preds3$upper2[2],
                                           preds3$lower2[6] - preds3$lower2[2])),
                                    0, #War =0
                                    1),  #officials = 1
                                  c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                           preds3$upper2[7] - preds3$upper2[3],
                                           preds3$lower2[7] - preds3$lower2[3])),
                                    1, #War =1
                                    0), #officials = 0
                                  c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                           preds3$lower2[8] - preds3$lower2[4],
                                           preds3$upper2[8] - preds3$upper2[4])),
                                    1, #War =1
                                    1) #officials = 1
                                  
))

colnames(po.rt.diff) <- c("LCI", "estimate", "UCI", "War", "PiS")

###
#TUSK
###
sim.df <- data.frame(conspir.tusk = sample(rep(c(0,1,0,1), length.out = 250)),
                     War = sample(rep(c(0,1,0,1), length.out = 250)),
                     officials = sample(rep(c(0,1,0,1), length.out = 250)),
                     log_follow = sample(seq(min(df$log_follow, na.rm = T), 
                                             max(df$log_follow, na.rm = T),
                                             length.out = 500)),
                     log_friends = sample(seq(min(df$log_friends, na.rm = T), 
                                              max(df$log_friends, na.rm = T),
                                              length.out = 500)),
                     verified = sample(c(rep(TRUE, 33), rep(FALSE, 467))),
                     month = c(rep(1:3,length.out = 42), 
                               rep(5:12, length.out = 112),
                               rep(4, length.out = 346)))

sim.df$conspir.tusk <- as.factor(sim.df$conspir.tusk)
sim.df$month <- as.factor(sim.df$month) 


preds2 <- cbind(sim.df, boot_pi(m.conspir.tusk.fav, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.tusk, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))

#Calculate differences
tusk.fav.diff <- as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                              preds3$upper2[5] - preds3$upper2[1],
                                              preds3$lower2[5] - preds3$lower2[1])), 
                                       0, #War =0
                                       0), #officials = 0
                                     c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                              preds3$upper2[6] - preds3$upper2[2],
                                              preds3$lower2[6] - preds3$lower2[2])),
                                       0, #War =0
                                       1),  #officials = 1
                                     c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                              preds3$upper2[7] - preds3$upper2[3],
                                              preds3$lower2[7] - preds3$lower2[3])),
                                       1, #War =1
                                       0), #officials = 0
                                     c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                              preds3$lower2[8] - preds3$lower2[4],
                                              preds3$upper2[8] - preds3$upper2[4])),
                                       1, #War =1
                                       1) #officials = 1
                                     
))

colnames(tusk.fav.diff) <- c("LCI", "estimate", "UCI", "War", "PiS")

#TUSK & Retweets

preds2 <- cbind(sim.df, boot_pi(m.conspir.tusk.rt, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.tusk, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))

#Calculate differences
tusk.rt.diff <-as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                            preds3$upper2[5] - preds3$upper2[1],
                                            preds3$lower2[5] - preds3$lower2[1])), 
                                     0, #War =0
                                     0), #officials = 0
                                   c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                            preds3$upper2[6] - preds3$upper2[2],
                                            preds3$lower2[6] - preds3$lower2[2])),
                                     0, #War =0
                                     1),  #officials = 1
                                   c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                            preds3$upper2[7] - preds3$upper2[3],
                                            preds3$lower2[7] - preds3$lower2[3])),
                                     1, #War =1
                                     0), #officials = 0
                                   c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                            preds3$lower2[8] - preds3$lower2[4],
                                            preds3$upper2[8] - preds3$upper2[4])),
                                     1, #War =1
                                     1) #officials = 1
                                   
))
colnames(tusk.rt.diff) <- c("LCI", "estimate", "UCI", "War", "PiS")


###
#Russia
#Note: In our original submission, we included analyses for the CT that Russia
#or Putin caused the plane to crash. We maintain the analysis below for full
#transparency and to fully replicate the simulation seeds for all subsequent models.
###
sim.df <- data.frame(conspir.rus = sample(rep(c(0,1,0,1), length.out = 250)),
                     War = sample(rep(c(0,1,0,1), length.out = 250)),
                     officials = sample(rep(c(0,1,0,1), length.out = 250)),
                     log_follow = sample(seq(min(df$log_follow, na.rm = T), 
                                             max(df$log_follow, na.rm = T),
                                             length.out = 500)),
                     log_friends = sample(seq(min(df$log_friends, na.rm = T), 
                                              max(df$log_friends, na.rm = T),
                                              length.out = 500)),
                     verified = sample(c(rep(TRUE, 33), rep(FALSE, 467))),
                     month = c(rep(1:3,length.out = 42), 
                               rep(5:12, length.out = 112),
                               rep(4, length.out = 346)))

sim.df$conspir.rus <- as.factor(sim.df$conspir.rus)
sim.df$month <- as.factor(sim.df$month) 


preds2 <- cbind(sim.df, boot_pi(m.conspir.rus.fav, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.rus, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))

#Calculate differences
rus.fav.diff <- as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                             preds3$upper2[5] - preds3$upper2[1],
                                             preds3$lower2[5] - preds3$lower2[1])), 
                                      0, #War =0
                                      0), #officials = 0
                                    c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                             preds3$upper2[6] - preds3$upper2[2],
                                             preds3$lower2[6] - preds3$lower2[2])),
                                      0, #War =0
                                      1),  #officials = 1
                                    c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                             preds3$upper2[7] - preds3$upper2[3],
                                             preds3$lower2[7] - preds3$lower2[3])),
                                      1, #War =1
                                      0), #officials = 0
                                    c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                             preds3$lower2[8] - preds3$lower2[4],
                                             preds3$upper2[8] - preds3$upper2[4])),
                                      1, #War =1
                                      1) #officials = 1
                                    
))

colnames(rus.fav.diff) <- c("LCI", "estimate", "UCI", "War", "PiS")


#RUS & Retweets

preds2 <- cbind(sim.df, boot_pi(m.conspir.rus.rt, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.rus, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))

#Calculate differences
rus.rt.diff <- as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                            preds3$upper2[5] - preds3$upper2[1],
                                            preds3$lower2[5] - preds3$lower2[1])), 
                                     0, #War =0
                                     0), #officials = 0
                                   c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                            preds3$upper2[6] - preds3$upper2[2],
                                            preds3$lower2[6] - preds3$lower2[2])),
                                     0, #War =0
                                     1),  #officials = 1
                                   c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                            preds3$upper2[7] - preds3$upper2[3],
                                            preds3$lower2[7] - preds3$lower2[3])),
                                     1, #War =1
                                     0), #officials = 0
                                   c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                            preds3$lower2[8] - preds3$lower2[4],
                                            preds3$upper2[8] - preds3$upper2[4])),
                                     1, #War =1
                                     1) #officials = 1
                                   
))
colnames(rus.rt.diff) <- c("LCI", "estimate", "UCI", "War", "PiS")


###
#Collusion
###
sim.df <- data.frame(conspir.collusion = sample(rep(c(0,1,0,1), length.out = 250)),
                     War = sample(rep(c(0,1,0,1), length.out = 250)),
                     officials = sample(rep(c(0,1,0,1), length.out = 250)),
                     log_follow = sample(seq(min(df$log_follow, na.rm = T), 
                                             max(df$log_follow, na.rm = T),
                                             length.out = 500)),
                     log_friends = sample(seq(min(df$log_friends, na.rm = T), 
                                              max(df$log_friends, na.rm = T),
                                              length.out = 500)),
                     verified = sample(c(rep(TRUE, 33), rep(FALSE, 467))),
                     month = c(rep(1:3,length.out = 42), 
                               rep(5:12, length.out = 112),
                               rep(4, length.out = 346)))

sim.df$conspir.collusion <- as.factor(sim.df$conspir.collusion)
sim.df$month <- as.factor(sim.df$month)  


preds2 <- cbind(sim.df, boot_pi(m.conspir.collusion.fav, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.collusion, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))

#Calculate differences
collusion.fav.diff <- as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                                   preds3$upper2[5] - preds3$upper2[1],
                                                   preds3$lower2[5] - preds3$lower2[1])), 
                                            0, #War =0
                                            0), #officials = 0
                                          c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                                   preds3$upper2[6] - preds3$upper2[2],
                                                   preds3$lower2[6] - preds3$lower2[2])),
                                            0, #War =0
                                            1),  #officials = 1
                                          c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                                   preds3$upper2[7] - preds3$upper2[3],
                                                   preds3$lower2[7] - preds3$lower2[3])),
                                            1, #War =1
                                            0), #officials = 0
                                          c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                                   preds3$lower2[8] - preds3$lower2[4],
                                                   preds3$upper2[8] - preds3$upper2[4])),
                                            1, #War =1
                                            1) #officials = 1
                                          
))

colnames(collusion.fav.diff) <- c("LCI", "estimate", "UCI",  "War", "PiS")


#COLLUSION & Retweets

preds2 <- cbind(sim.df, boot_pi(m.conspir.collusion.rt, sim.df, 1000, 0.95))
preds3 <- as.data.frame(preds2 %>% group_by(conspir.collusion, War, officials) %>%
                          summarize(across(pred:upper, mean)))
preds3$pred2 <- exp(log(preds3$pred))
preds3$upper2 <- exp(log(preds3$upper))
preds3$lower2 <- exp(log(preds3$lower))

#Calculate differences
collusion.rt.diff <- as.data.frame(rbind(c(sort(c(preds3$pred2[5] - preds3$pred2[1],
                                                  preds3$upper2[5] - preds3$upper2[1],
                                                  preds3$lower2[5] - preds3$lower2[1])), 
                                           0, #War =0
                                           0), #officials = 0
                                         c(sort(c(preds3$pred2[6] - preds3$pred2[2],
                                                  preds3$upper2[6] - preds3$upper2[2],
                                                  preds3$lower2[6] - preds3$lower2[2])),
                                           0, #War =0
                                           1),  #officials = 1
                                         c(sort(c(preds3$pred2[7] - preds3$pred2[3],
                                                  preds3$upper2[7] - preds3$upper2[3],
                                                  preds3$lower2[7] - preds3$lower2[3])),
                                           1, #War =1
                                           0), #officials = 0
                                         c(sort(c(preds3$pred2[8] - preds3$pred2[4],
                                                  preds3$lower2[8] - preds3$lower2[4],
                                                  preds3$upper2[8] - preds3$upper2[4])),
                                           1, #War =1
                                           1) #officials = 1
                                         
))
colnames(collusion.rt.diff) <- c("LCI", "estimate", "UCI", "War", "PiS")


comb <- rbind(po.fav.diff, po.rt.diff, 
              tusk.fav.diff, tusk.rt.diff,
              collusion.fav.diff, collusion.rt.diff)

comb$model <- c("PO", "PO", "PO", "PO", "PO", "PO", "PO", "PO", 
                "Tusk", "Tusk", "Tusk", "Tusk",  "Tusk", "Tusk", "Tusk", "Tusk",
                "Collusion", "Collusion",  "Collusion", "Collusion", "Collusion", "Collusion",  "Collusion", "Collusion")
comb$War <- as.factor(comb$War)
comb$engagement <- c("Likes", "Likes","Likes","Likes",
                     "Retweets", "Retweets", "Retweets", "Retweets", 
                     "Likes", "Likes","Likes","Likes",
                     "Retweets", "Retweets", "Retweets", "Retweets", 
                     "Likes", "Likes","Likes","Likes",
                     "Retweets", "Retweets", "Retweets", "Retweets")
comb$War2 <- ifelse(comb$War == 1, "After February 24", "Before February 24")
comb$PiS2 <- ifelse(comb$PiS == 1, "Yes", "No")
comb$PiS3 <- ifelse(comb$PiS == 1, "PiS Official", "Non-PiS Official")



###
pd <- position_dodge(width = 0.6)

trip.inter.plot <- ggplot(comb, aes(x = engagement, y = estimate, shape = factor(PiS2), color = War2)) + 
  geom_point(position = pd, size = 1.5) + 
  coord_flip() +
  scale_color_grey(start = .2, end = .5, name = "War?") +
  geom_errorbar(aes(ymin = LCI, ymax = UCI), width=.1, position = pd) +
  labs(y = "Difference in the Predicted Level \n of Engagement When Content Included",
       x = " ",
       shape = "PiS Official?")+  
  facet_wrap(~model, ncol = 1) +
  ylim(-180, 75)  +
  
  theme_classic(base_size=20) +
  #theme(legend.position = "bottom") +
  geom_hline(aes(yintercept = 0), color = "blue", linetype = "dotted")
trip.inter.plot
